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Abstract This paper considers mutual obligations in the interconnected bank 
system and analyzes their influence on joint and marginal survival probabilities as 
well as CDS and FTD prices for the individual banks. To make the role of mutual 
obligations more transparent, a simple structural default model with banks’ assets 
driven by correlated multidimensional Brownian motion with drift is considered. 
This model enables a closed form representation for many quantities of interest, 
at least in a 2D case, to be obtained, and moreover, model calibration is provided. 
Finally, we demonstrate that mutual obligations have to be taken into account in 
order to get correct values for model parameters. 

Keywords 2D structural default model, mutual obligations, joint and marginal 
survival probabilities, CDS and FTD prices 


1 Introduction 

Structural default framework is widely used for assessing credit risk of a corporate 
debt. It was introduced in a simple form in the seminal work [21], and was further 
extended in various papers, see a survey in [20] and references therein. In con¬ 
trast to reduced-form models (see, e.g., [5]) structural default models suffer from 
the curse of dimensionality when the number of counterparties grows; however, 
these models provide a more detailed and financially meaningful description of the 
default event for a typical firm. 
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Inspired by [25], in [15] we extended the structural framework by taking into 
account the fact that banks have mutual liabilities among themselves. Accounting 
for this effect is very important in order to accurately analyze credit worthiness of 
individual banks and the banking network as a whole. For instance, large mutual 
liabilities imply that adverse shock to a bank is rapidly transmitted to the entire 
system, with severe implications for its stability ([8]). The authors of [8] indicate 
that renegotiations between highly interconnected banks facilitate mutual private 
sector bailouts to lower the need for government bailouts. The relative size of 
mutual liabilities compared to total liabilities is fairly substantial. For instance, 
the relative fraction of interbank loans is 12% in the EU, 8.5% for Canada ([8]), 
and 4.5% for US (as per Economic Research website of the Federal Reserve Bank 
of St. Louis). 

An extended Merton model with mutual liabilities and continuous default mon¬ 
itoring can be built by combining correlated Merton balance sheet models, cali¬ 
brated by using observed bank equity returns, and analyzing potential clearing of 
interbank liabilities in the spirit of [9], In [15] we assumed that banks’ assets are 
driven by correlated Levy processes with idiosyncratic and common components 
and developed a novel pseudo partial differential equation computational method 
in order to make the problem of computing joint and marginal survival prob¬ 
abilities tractable. The effect of mutual liabilities was discussed, and numerical 
examples were given to illustrate its importance. 

Obviously, the knowledge of joint and marginal survival probabilities is impor¬ 
tant for successful calibration of the model to Credit Default Swap (CDS) spreads 
and first-to-default (FTD) instruments. Since the general case is fairly compli¬ 
cated, here we restrict ourselves to the case when banks’ assets are driven by 
correlated Brownian Motions with drifts. Then in the 2D case we obtain explicit 
expressions for several quantities of interest including joint and marginal survival 
probabilities as well as CDS and FTD prices. Despite the fact that the model 
under consideration does not incorporate jumps, it is still beneficial as it enables 
an analytical assessment and provides a natural link to the analytical framework 
considered in [19, 20]. 

The rest of the paper is organized as follows. In Section 2 we propose a model 
for the general case of N banks. Sections 3,4 present the governing equations and 
Green’s function approach to the solution of these equations for joint and marginal 
survival probabilities for two banks with mutual obligations. In Sections 5,6 the 
prices of CDS and FTD contracts are calculated, and results of our numerical 
experiments are presented. We also validate our results by comparing analytical 
solutions with numerical solutions obtained by using a finite difference algorithm 
described in [15]. Section 7 discusses calibration of the model presents some nu¬ 
merical results. Our conclusions are presented in Section 8. 


2 Model 


Consider a set of N banks with external assets and liabilities A,, L,;, i = 1, ...,7V, 
and interbank assets and liabilities Lji, j = 1,..., N, respectively. In other words, 
Lij is the amount the i-tli bank owes from the j-th bank, etc. Thus, total assets, 
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liabilities and capital of the i-th bank have the form 


-'1 / — -1 i “ 1 “ ^ ) Lji • Li — Li + ^ ' Li 
jAi 


Ci — 1 -/ Li — -1 i A i 


where 

Ai = Lj + y, (Ljj — Ljj). (1) 

i^i 

For simplicity, we assume that the corresponding dynamics is governed by the 
SDEs of the form 


dA it — rA it dt + , (2) 

dL.i t rLi t dt, dLijj i Lijjdt, 

subject to the initial conditions A^ 0 = A,(0), Li j0 = Li(0), Lij y0 = (0), so that 

L-i f , Lij t are deterministic functions of the time t 1 . In Eq.(2) r is the risk-free rate, 
<ji is the volatility of the i-th asset (which is assumed to be constant), and W^t 
is the corresponding Brownian Motion. Elements of the corresponding correlation 
matrix are denoted by pij. The above assumptions can be generalized in a variety 
of ways, which will be discussed elsewhere. 

We assume that all the liabilities (both external and interbank) are settled at 
a certain maturity T > 0. Thus at t = T, the i-th bank defaults if < Lj(T), 
or, equivalently, if A, ( t < Xf(T). Below we denote the default time of the ith bank 
by Tj. 

The k-th bank defaults at < T. We describe defaults at intermediate times 0 < 
Ti < T in the spirit of [6] by assuming that the i-th bank defaults at time Ti 
provided that 


A i(n) < (3) 

A i (Ti) — Hi ^T ? ;(Tj) T ^ Ljj (T,)j ^ Lji(Ti) : 

where 0 < Rj < 1 is the recovery rate, which is assumed to be constant up to the 
time t = T. We need to emphasize that according to these settings, the default 
boundary is discontinuous at t = T, because Rj experiences a jump at this point 
from its value Rj at t < T to 1 at t = T (and so \f transforms to A=). 

These default boundaries are valid if no other bank defaults till t = T. Now 

assume the opposite, i.e. that the fc-th bank is the first to default at time r k < T, 

and so we are left with a set of N — 1 surviving banks. At time the assets and 
liabilities of the i-th bank, i k, have the form 

Ai(Tj~) Ai(r k ) ^ Lji(r k ) + RkL k i(r k ), 

Li(T k ) = Tj(rfc) + ^ Lij(r k ) Lj k (T k ). 


1 Accordingly, further we will denote them as Li ( t ), Lij (t). 
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We assume that for surviving banks mutual liabilities stay the same, while 
their external liabilities jump according to the rule 


Lii^k) ^ Lj J (rk') — Ti{l~k) T ^airk) RfcL k i i^k )■ 

Surviving banks’ capital naturally takes a hit 

Ci{rk) ^ Ci( t ^) = Ci^Tk) (1 Rk)Lki{j~k)- 

Thus, each default reduces the set of surviving banks and modifies the corre¬ 
sponding default boundaries as 


Aifc(t) — 


_ I t <T, 


i k 


t — T, 

A^(t) = Ri[Li(t) + Lik{t) — RuLki(t)\ + [RiLij(t) — Lji(t)\. 

Kk(T) = Li(t) + Lik(t) — RkLki(t) + [Lij(t) — Lji(t )], 


( 4 ) 


As an example, consider N = 2. Then for the remaining bank k = 3 — k, the 
default boundary is given by ([15]) 2 


f R k (-kfc + L k k ~ RkL kk )» Tk<t <T, 
\L k + L k k ~ Rk^kk’ T k <t = T. 


It is clear that 


k ~ ~ “ 


f (1 — R k Rk) L kk = , 

1(1 ~Rk)L k - k = A\^ 7 


t = T, 


( 5 ) 


and AXj, > 0. 

To make these definitions more transparent the computational domain is repre¬ 
sented in Fig. 1. Here, if there are no defaults, we have a rectangular computational 
domain which lies above the piece-wise constant line 5 — 3 — 4. If the bank 2 de¬ 
faults, this domain transforms to that which lies to the right of the line 5 —3 —7 —8. 
If the bank 1 defaults, the domain transforms to that which lies above the line 
1 - 2 - 3 - 4. 


The i-th bank defaults at Ti = T. In this case the definition of A i in Eq.(5) should be 
changed. Indeed, if assets of the i-th bank breach below its liabilities at some time 
before maturity, the bank has some period of time to recover, unless it breaches 
below the level \f . At this level the bank’s counterparties don’t believe anymore 
in its ability to recover, and it defaults. Obviously, at t = T the bank doesn’t have 
time to recover. Therefore, at the most it can pay to its obligors the current amount 
of money in hands, i.e. the total value of the bank assets 3 which is a fraction 
7i, 0 < 7j < 1, of its liabilities. Accordingly, in the spirit of [9], to determine 

2 In this case we omit the second index of A,/, as the defaulted bank is determined uniquely. 

2 We consider an idealistic situation when all bank’s assets upon default can be immediately 
converted to cash with no delay and further losses. 
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Fig. 1: Default boundaries of two banks with and without mutual liabilities at t < T. 


default boundaries we need to find the vector 7 = { 7 j}, 0 < 7 , < 1, i € [1, JV] 
which solves the following piece-wise linear problem in the unit cube: 


min \ Ai(T) + '£ / 'y J L ji (T), L^T) + £ L tJ (T) \ = 7i ( L t (T) + Y L n( T ) . 








( 6 ) 

Introducing new non-dimensional variables = Ai(T)/ Li(T), lj t = Lji[T)/Li(T) 
the problem given in Eq.( 6 ) can be re-written in the form 


min < a, + Ijljii 1 / = 7 i- (7) 

{ i& ) 

It is clear that 7 ; = 1 (so that the i-th bank survives) if ai + > 1. 

And 7 j < 1 otherwise, so that the i-th bank defaults. This description suggests 
that defaults in the interlinked set of banks can happen outright, when 

Ai ( T) < X= (T), 


and through contagion, when 

Li(T) + Y, [L i:j (T) - 7 jLji (T)] > Ai (T) > X= (T). 












6 


Andrey Itkin, Alexander Lipton 


Eq.(7) can be uniquely solved. A brief discussion is given in Appendix A. 
Accordingly, in this case we change the definition of (T) at 77 = T to 

Kt = Li(T) + J2 Ly (T) - Y, R j T (^)Lji(T) (8) 

*7 H 


where 


Rj,r( 7 ) = min 1 , 


_ ijH _ 

Lj(T) + Y^ L R(T) 


7 = [71,-,7 a ]- (9) 



Fig. 2: Default boundaries of two banks with and without mutual liabilities at t = T. The dot 
pattern marks the whole computational domain X). 


It follows that the default boundary \,,t piece-wise linearly depends on all 
Aj(T), j £ [1,1V], j i. In particular, let IV = 2 and T 2 = T, hence when 
A 2 (T) = 0 we have from Eq.( 8 ) 

Ai t = 7—.—7—, A = L1L2 + L12L2 + L\Ij 2 i. 

T 2 + A 21 


Ai — Ai t = L2i{R.2 t(1) — R 2 ) = L 21 

r 2 <T 


L 2 + L 2 



Therefore. 
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This behavior is illustrated in Fig. 2. Since tj = T, and, thus, Rj = 1, from 
Eq.(5) we have A= = A=, or AXj = 0. Therefore, when ^(T) grows from 0 to A^, 

the default boundary A^t moves from ^ + L 21 ^ along the line 8-9-6. 

At point 6 the default boundary Xi t transforms to A^ = Af, and further 
doesn’t depend on A^(T) when the latter increases. This occurs at the point 
A 2 = (T) = A^. Thus, the whole default boundary of the first bank in Fig.2 can 
be seen as a line passing through the points 8-9-6-5. Similarly, for the second bank 
the default boundary in Fig. 2 can be seen as a line passing through the points 
1-2-6-4. 

Also in Fig. 2 D 12 is the domain where both banks don’t default, Di - where 
the second bank defaults while the first one does not, D 2 - where the first bank 
defaults while the first one does not, D is the whole computational domain marked 
by the dot pattern, and in the domain D = S)\Di 2 both banks default. 

As always, it is useful to describe the evolution of the set of banks under 
consideration in terms of non-dimensional variables. To this end, we introduce the 

average volatility c o = a i) , and define 


t = u t, T = uT, 


X if = - hr 


t _ 

s i — o 

2 to 


The corresponding dynamics of X i j is governed by the SDE: 

d Xi,i = iidt + dWij, 

while the default conditions now transform to X. t < fxi, with fj,i defined as 


f Hi = 0, 

^ ® = \ u = = M. l n (K® 


By definition, >0. 

Below we omit bars for the sake of simplicity. 

If the j- th bank defaults at Tj < T, then for the i- th bank the default boundary 
is given by 




i< = M In l 


= _ u> (\jip) 


r,= = At ln XI 

~ tJj 111 1 x< 


a mj’ 


Note, that according to Eq.(2) /lA doesn’t depend on t. 

It can be seen that the boundary condition in Eq.(12) at t = T doesn’t match 
to the terminal condition which, according to Eq.(8), reads 


= — ln 


A ij,T( T ) \ , 

A fj{T) ) ^ 


Mathematically, this means that our problem belongs to the class of problems 
with a boundary (transition) layer at t = T. Financially, the behavior of the 
solution in this layer is determined by the detailed specification of the contract. 
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For instance, if the bank is close to maturity, say a day before, the recovery rate 
could be defined to smoothly transit from f?j to 1 within this last day. Or, some 
other conditions specific to the contract in question could be issued. However, 
we don’t consider these details, assuming that the boundary layer is thin, and, 
therefore, any perturbation of the solution due to the existence of this layer, is 
dumped out pretty fast when moving away from this layer. In other words, as we 
ignore a detailed consideration of the boundary layer, our solution experiences a 
jump at t = T. Therefore, after the closed form solution is found we will compare 
it with the numerical solution of this problem to reveal sensitivity of the former 
to the value of the described effect. 

Below we provide all the results just for two-dimensional case N = 2 while the 
multi-dimensional case will be presented elsewhere. Accordingly, in the definition 
of jj.fj for easiness of reading we will omit the second index as in this case it 
doesn’t bring any confusion. 


3 Governing equations 

Based on the analysis presented in the previous section, the joint survival prob¬ 
ability Q(t, Xi,X 2 ) of two assets Xi,X 2 is defined in the domain i?(t,Xi,X 2 ) : 
[0,T] x [0, 00 ] x [0, oo] 4 5 . It solves the following terminal boundary value problem 
([ 20 ]) 


Qt ( t,Xi,X 2 )+CQ (t, Xi,X 2 ) 
Q (T, A'i, A 2 ) = lygDij, Q(t, 0,A 2 ) 


0, 

0, Q (t, Ai, 0) = 0, 


where 


(14) 


CQ — 4\pQ + £ • VQ, 

_ 1 d 2 d 2 Id 2 

p = 2dXf +P dX 1 dX 2 + 2<9Xj’ 


£ — (£i,£2) t , 


l x is the Heaviside step function defined with the half-maximum convention , 
and the area Di 2 is defined in Fig. 2. We emphasize that the domain D; in X 
variables has a curvilinear boundary which depends on the value of Aj(T). Indeed, 
based on the definitions in Eq.(3), Eq.(8), Eq.(9) and Eq.(13), one can find, e.g., 
for i = 1 


Mi ,T 


= — In 
<ri 


L 1 + L\2 — R2,t( 1)T 2 i 
L 1 + L \2 — L21 


— In 
ci 


1 + 


T 2 i(l — -R 2 ,t(1)) 
L 1 + L 12 — L 21 


Next we define the corresponding marginal survival probabilities qi (t, Xi,X 2 ), 
i = 1,2, which are functions of both X\ and X 2 , also in the domain f2(t,X i,X 2 ). 
For brevity we provide all definitions and formulae for the first bank ( i = 1) while 

4 The space sub-domain of 12 corresponds to the dotted area in Fig. 2. 

5 Since the detailed consideration of the transition layer at t = T is omitted, this condition 
allows getting the correct value of xi(T p-T)’ see Eq.(22). 
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for the second one it could be done by analogy. So q\ (t, X \, A' 2 ) solves the following 
terminal boundary value problem 


qi,t ( t , X U X 2 ) + Cq 1 (t, X lt X 2 ) = 0, 
qi (T,X 1 ,X 2 ) = 1 xc. i), 2 un, ji 

qi (t, 0, X 2 ) = 0, Q1 (t,X i,0) = S(t,X 1 ) = 


fxi,o(f, A'i), 

\0, 


(15) 


Xi > sf, 

0< Ai <Ai, ’ 


9i (t, A'i, A 2 t 00 ) — Xi,oo(b A'i), gi (t, A'i t 00 , X 2 ) — 1. 


In Eq.(15) the domain Di is defined in Fig. 2. Function xi,o(t,Xi) is the ID 
survival probability, which solves the following terminal boundary value problem 


dtxi,o(t, x i) + £ixi,o(t, x i) = 0, (16) 

Xi,o(T,Xi) = l Xl>p =, Xi,o(f,Af) = 0, 

where 

C . = ±#_ + C.J_ 

1 2 i)Xf 

Accordingly, function xi,oo(t,X\) is the ID survival probability, which solves 
the following terminal boundary value problem 


dtXi,oo(t,Xi) + £iXi,oo(f,A'i) — 0, (17) 

Xi,oo(7\ X\) = l Xl>p =, Xi,oo(^,0)=0, 


4 Survival probabilities 

We solve Eq.(14) and Eq.(15) by introducing the Green’s function G(t, A'i, X 2 |t / , X \, 
A' 2 ), where A'(, X' 2 are the initial values of A'i, X 2 at t = t'. Below, where it is not 
confusing, for brevity we will also use the notation G(t — t', Xi,X 2 ), thus explicitly 
exploiting the fact that for our problem the Green’s function depends only on 
t — t', and omitting the second pair of arguments. The Green’s function solves the 
following initial boundary value problem 


Gt(t-t',X 1 ,X 2 )-tfG(t-t',X 1 ,X 2 ) = 0, (18) 

G( 0, A'i, X 2 ) = S(X 1 - X[)5(X 2 - X 2 ), 

G(t-t', 0,A' 2 ) = 0, G(t — t', A'i, 0) =0, 


where C) = A p — £ • V. A simple calculation yields 

(QG) t + CGQ - QC^G = 0, 

or, explicitly, 

(QG), + v(l iQx - G - QGx ‘ > - ^ + bQG ) = 0 

\l(Qx,G-QGx,) + liQx,G + (2QG) 
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(19) 


( 20 ) 


The Green’s theorem ([17]) yields 

rOO pOO 

Q(t',X[,X ' 2 ) = dX i/ dX 2 G(T-t',Xi,X 2 ) 

Jo Jo 

IJ G (t',X 1 ,X 2 ) dXidX 2 , 

(x 1 ,x 2 )eD 12 

where 7 =T — t'. Similarly, 

qi(t ', X\ ,X 2 ) = IJ G(t',X i, X 2 ) dX i dX 2 

(X 1 ,X 2 )e[D 12 UD 1 ] 

, OO 

+ 11' ds j dX 1 Gx 2 (r , - S ,X 1 ,0) X i,o(s,X 1 ) 

OO 

- J dXiGx 2 (t - s, X 1,0) XI, OO (s,Xi) 
o 

We start with noting that the ID Green’s function gi(d, Xi), d = t — t' at 
X 2 (t) < p 2 has the form 


gi(tl, XT) = 
Accordingly, 

Xi,o(t',X[) = / e“ 


s -^tf/2+| 1 (x 1 -x() 

77rd 


(X 1 -X [) 2 (x 1 +x[-2fi<Y 

g 2i9 — g 


( 21 ) 


OO 

(Xj-Xj) 2 

( xi+ x ;-2ii<) 2 ■ 

J e -e i r'/2+Zi(X 1 -X , 1 ) 

e 2 t' 

e 27 

V2ttt' 

V2ttt' 

hr 

- 

- 


dXi (22) 


-/ 

Ar 

= X | 
and 


oo (Xl -x; -f 1T ') 2 

e J7 7 


Tfci 


dXi -e 


OO ( j 2 y ^) 2 

- 2 { 1 (X(-/i<) f e 


77 


dX i 


Mr 


77 


f Ar + X[ - 2p< - Cl/ 


X - 


Xi,oo(*',xO = iV( - 


77 


hr - - 77 \ _ e - 2€l xj JV / hr + x( - Ci7 


^ t. _/ 

77 


77 


The corresponding 2D Green’s function has the form (see [18, 26] and references 
therein) 


G(0,Xi,X 2 ) = 

2 _(£77 +< x-A-', e) -i^+^ 2 


w&p 


rr! 


2d ^ ii/„ I —r- ) sin ( v n (j >) sin (i/„<7) 


(23) 
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where (,} denotes the dot product, Ik(x) is the modified Bessel function of the 
first kind, 


C = 


vo = < 


1 P 

P 1 


C -1 = J_ ( 1 ~P 


~P 1 , 


Q _ /~1~ 1 £ _ WK 

6 — C £, V n — , 

VO 


[" 7T + arctan (—p/p ), p > 0 
tt/2, 


p = 0, R 2 = {X , C~ 1 X T ), R' 2 = 


y arctan (— p/p ), p < 0, 

' ir + arctan ^ x , 2 ) > A 2 < pAd, 

<5(Ai, X 2 ) = { 7r/2, X 2 = pXi, 


arctan 


pAi \ 

—pA'i + X 2 ) ’ 


A 2 >pAi, 


0 = $(X 1 ,X 2 ), <t> = $(X'i,X' 2 ), X = (Ai, X 2 ). 


-2 -i 2 

p = i - p , 
(x 7 , C~ 1 X ,T ), 


Accordingly, 






. J2 (-1)” +1 Union sin (" 0 


(24) 


Substitution of these formulas into Eq.(19), Eq.(20) yields senri-analytical ex¬ 
pressions for Q and qi. However, from the computational point of view, it is more 
efficient to introduce a new function 


qi(t, Xi,X 2 ) — qi(t,X i,X 2 ) — Xi,oo(t, Ai). 

In contrast to qi(t, Xi,X 2 ) this new function solves a problem similar the problem 
given in Eq.(15), but with a homogeneous upper boundary condition: 

qi,t(t, Xi, X 2 ) + Cqi(t, Xi, X 2 ): = 0, (25) 

<h(t,0, A 2 ) = gi(t, Ai,X 2 f oo) = 0, 

9i (T, X'i, X 2 ) = — ljcgD> 

where D is the area inside the curvilinear triangle with the vertexes in points 6-7-9 
in Fig. 2. 

As the equations in Eq.(15) and Eq.(25) differ just by the source term, the 
Green’s function of Eq.(25) is also given by Eq.(23). Accordingly, the solution of 
Eq.(25) reads 

q 1 (t , ,X' 1 ,x! 2 )=xi,oo(t^X' 1 )- jj G(r , ,X 1 ,X 2 )dX 1 dX 2 (26) 

+ lf Q dsj^ Gjr a (T , -a,X 1 ,0)[S(T , -«,X 1 )-x ll oo(T , -s,X 1 )]dX 1 . 

Another simplification could be made if one wishes to compute the first in¬ 
tegral in Eq.(20). For a better accuracy it could be reasonable to represent it as 
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a difference of two integrals. The first one is taken over the positive quadrant 
(Xi,X 2 ) £ [0, oo) x [0, oo) while the second integral is defined in a union of two 
semi-infinite strips: T>i U T >2 U V. The idea is that the second integral is defined 
in the area which is finite either in one or the other direction, and the first inte¬ 
gral could be represented in the closed form. Therefore, the total computational 
error is less. We underline, that to the best of our knowledge representation of the 
first integral in closed form yet was not given in the literature, so we present this 
derivation in Appendix B. 


4.1 Numerical experiments 

In our test examples we solved Eq.(15) by using first a finite difference scheme (FD) 
and then compared it with the analytical solution given by Eq.(26). Since Eq.(15) 
is a pure convection-diffusion two-dimensional problem, we solved it numerically 
by using a Hundsdorfer-Verwer scheme, see [12]. A non-uniform finite-difference 
grid was constructed similar to [14] with the grid nodes concentrated close to 
/},= , * = 1,2. We solved the problem using parameters given in Table l 6 : 


Table 1: Parameters of the structural default model. 


L i,o 

^2,0 

hl2,0 

^21,0 

Ri 

R 2 

T 

^1 

O- 2 

p 

60 

70 

10 

15 

0.4 

0.45 

1 

1 

1 

0.5 


We computed all tests using a 100 x 100 spatial grid. Also we used a constant 
step in time At = 0.01, so that the total number of time steps for a given maturity 
is 100. The marginal survival probability X 2 ) at t = 0 computed by using 

this method is presented in Fig. 3. 



Fig. 3: The marginal survival proba¬ 
bility qi(Xi,X 2 ) computed by using a 
Hundsdorfer-Verwer scheme. 


Fig. 4: The difference between marginal 
survival probabilities qi(Xi, X 2 ) computed 
with and without mutual obligations using 
the FD method. 


6 In our setting the value of the interest rate r doesn’t matter. 
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It is easy to see that for our chosen parameters jLf = 0.6659, = 0.2548, /q = 

1.4424, /if = 0.9764, flf = 1.5821, /if = 1.0534. 

To observe the effect of the mutual liabilities we repeated this test, but with 
zero mutual liabilities. Therefore, as compared with the previous case, now the 
total assets of the i-th bank are Ai + JT while its liabilities are Li + JT Lji. 
But to provide a correct comparison we need to keep the asset values Ai constant. 
Therefore, in this case we re-adjust liabilities to Li + JT Lji — JT . In words, 
that means that, if JT Ljj is positive, the bank i gets extra cash and then spends 
it retiring some of its external liabilities. If this amount is negative, then it is 
borrowed from the external sources. After this adjustment is done we set L 21 = 
L 12 = 0 in our calculations. In what follows we call this procedure an Adjustment 
Procedure (AP). 

The difference of two solutions is presented in Fig. 4. The above plot clearly 
demonstrates a significant difference in the solution in the area close to AT = 
/if, A '2 = /rf, i.e. the effect of the mutual liabilities is pronounced in this area. 

Next we want to compare the analytical and FD solutions. Since the integrands 
in Eq.(26) are highly oscillating functions, to get a reasonable accuracy we used 
a Gauss-Kronrod algorithm in both directions. Fig. 5 demonstrates the difference 
in these solutions when there are no mutual obligations. Both solutions coincide 
pretty well. However, when mutual obligations are taken into account the difference 
increases as it can be seen in Fig. 6. The difference is bigger in the area closer to 
A] = /if ■ 



Fig. 5: The difference between marginal sur¬ 
vival probabilities qi(Xi, X 2 ) computed by 
the analytical and FD methods with no mu¬ 
tual obligations. 


Fig. 6: The difference between marginal sur¬ 
vival probabilities qi(X\, X 2 ) computed by 
the analytical and FD methods with mutual 
obligations, T = 5 years. 


Performance-wise, computation of the marginal probabilities on this spatial 
grid using the FD scheme takes 21 secs for T = 1 year. At the same time, compu¬ 
tation of a single point on a grid using our analytical methods takes 0.4-0.6 secs 7 . 
Therefore, if the marginal probabilities should be computed at every node of the 
FD grid, using the analytical method it would take about 4000 secs, which is pretty 
slow. However, when calibrating the model with unknown volatilities and the cor¬ 
relation coefficient, we need functions at only three quotes that could be taken 

7 We ran this test in Matlab on a standard PC with Intel Xeon E5620 2.4 Ghz CPU. 
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from the market values of the CDS spreads and First-To-Default swaps spreads. 
Thus, in our simplistic model we need just three points, which takes about 1.2 secs 
to compute using our analytics. In contrast, the FD scheme cannot be reduced just 
to three points on the grid, and, therefore, for such kind of calibration is much less 
efficient than the analytical method. This is the reason we propose the approach 
of this paper for doing fast calibration of the model. 

In a more general setting, e.g., the one proposed in [15] this simplified approach 
could be used to produce a ’’smart” initial guess for the parameters of the marginal 
distributions. Then, using this guess, the whole rather complicated problem could 
be calibrated much faster than starting with some arbitrary values of the param¬ 
eters, since in this case only relatively small increments of the initial guess should 
be found. 


5 Pricing CDS contracts 

We now describe how to price CDSs and FTD in our setting. 


5.1 CDS 


The price of a CDS Ci(t, Xi, X 2 ) 8 written on the first bank solves the following 
problem, ([5]): 


<7 M (t,Xi,X 2 ) + CCi(t,X 1 ,X 2 ) = ft, 


Ci(t,X 1 ,0) = ¥(t,X 1 ) 


fci,o(t,Xi), X\ > jj,f 

\l-Ri, Xi </*<’ 


<7i(t,0,X 2 ) — 1 — Ri, C 1 (t,X 1 ,X 2 f 00 ) — ci,oo(t,Xi), 


(27) 


where c* is the coupon rate, ci,o(t, AT) is the solution of the one-dimensional ter¬ 
minal boundary value problem 


dtc\ t o(t, AT) + C\c\ t o(t, Xi) — ?i, ( 28 ) 

ci,o(t, fri ) = 1 - Ri, ci,o(t, 00 ) = — ?i (T - t), 

Cl,o(T,Xi) = (1 - 

and ci >0 o (t, XT) is the solution of another one-dimensional terminal boundary value 
problem 


dtCi t oo(t, Xi) + C\c\ t oo(t, Xi) = ci, (29) 

ci,oo(t,0) = 1 - Ri, ci,oo(t,oo) = -?i(T - t), 

ci,oo(T,XT) = (1 — Ri)lx 1 </ij= ■ 

Also the statement of problem given in Eq.(27) must be supplied with the 
terminal condition Ci(T,XT,X 2 ), which could be provided based on the picture 


For C 2 {t, X±, X 2 ) similar expressions could be provided by analogy. 
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presented in Fig. 2. Omitting some intermediate algebra, we obtain the following 
condition 


Ci(T, X \, X-2 ) — a il(x 1 ,x 2 )e[DuD 2 ] 


fl - min[^ 1)T (l), Rx], 
- min[7? 1)T (7 2 ),i?i], 


(Xi,x 2 ) e d 2 
(x u x 2 ) e D. 


The value of the components 7 , at (Xi,X 2 ) £ D are determined by solving the 
detailed balance equations which follow from the general IV-dimensional problem 
given in Eq.( 6 ) 


Ai + 7 2 L 2 i — 71 (L1 + L12 ), 
A 2 + 7l 1/12 = 72 (T 2 + I/ 2 l) , 
The solution in the original variables reads 

Lj Ai + LjJAi + Aj) 


7 i = 


A 


i = 1 , 2 . 


Observe that the Green’s function for Eq.(29) is that given by Eq.(21). There¬ 
fore, 


ci,o(t', X[) = (1 - Ri) r <? 1 (r / , Xi)dXi + ^ — [ 
Jjxf 2 Jo 


1-Ri f T dgi(r' -s,Xi) 


dXx 


ds 


Xi =Ai 


- ?i 


n oo 

9i(t - a, X^dX^s = h + I 2 + I 3 . 

;< 


(30) 

/o •/£< 

All these integrals can be computed in the closed form. Omitting some inter¬ 
mediate algebra we provide just the final results: 

h = (1 - Ri){ [N (-y-) - N (—2y_ -z)] + N (y+) - N (z) 


h = (l~Ri) [e- 2 ?l( - Y ^ A i <) TV( 2 /_) + TV(-y + )] , 


(31) 


h = SIT 


1 - -k±=N(-y+) - 


-2{ 1 (.Y 1 '-ii 1 < ) y- 


_ ±(A"i — yi ) + Cl t _ X'i — yi + t£ 1 

?/± - /- 5 ^ - 


Vt' 


W 1 


By analogy 


fMi 


ci,oo(t',-Ai) = (1 - 7?i) / ^(r'.XiJdXi + 


I-/R 1 f T dgi(r' — s,Xi) 
2 In dX 1 


ds 


x 1 =o 


Cl 


n oo 

gi(r' — s,X\)dX\ds, 


(32) 


where gi(r / ,AT) can be obtained from gi(r',Xi) by setting in Eq.(21) yf = 0. 
Accordingly, these integrals in closed form are given by Eq.(31) by replacing /xf = 0 
and yff = yf. 
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Using the same trick as in the previous section when we computed the marginal 
survival probability qi(t,Xi,X 2 ), the final solution of this problem could be rep¬ 
resented as follows: 

roo /*oo 

Ci(jf,X' 1 ,X' 2 ) = ci, 00 (t',X' 1 )+ / <t>(X 1 ,X2)G(T , ,X 1 ,X 2 )dX 1 dX2 (33) 

Jo Jo 

+ \ f ds [ Gx 2 (t — s,Xi, 0) \F(t — s, A'i) — ci l oo(r / — s, X\)\ dX\, 
z Jo Jo 

<t>{X i,X 2 ) = a i 1 (j ( : li A' 2 )e[DuD 2 ] — (-*- — 

where the Green’s function G (t 1 ,Xi,X 2 ) is given in Eq.(23). 


5.2 Numerical experiments 

It is well-known in a theory of heat conduction that a direct implementation of 
Eq.(33) is still impractical. The reason is that at X' 2 = 0 the first integral in Eq.(33) 
vanishes, so the second one must converge to ci j o(t , 1 X[) — ci i00 (t', X[) to provide 
the correct boundary condition at X 2 = 0. However, as it could be checked, at 
X 2 = 0 we have Gx 2 (t 1 — s,X 1 ,0) = 0, and, hence, the formal validation of Eq.(20) 
fails at the boundary. It is explained, e.g. in [16], the reason is that the series in the 
second integral in Eq.(33) is not uniformly convergent at X 2 = 0, so the transition 
to the limit X 2 —> 0 using this representation is complicated and impractical from 
the computational point of view. 

This problem, however, can be overcome by applying another elegant trick that 
we describe in more detail in Appendix C. 

Further on, we ran the same test as in the previous section with parameters of 
the model given in Table 1, and = 0.02. We used the same FD method to verify 
our solution, see the previous section for the description of the method. The CDS 
prices for t = 0 are presented in Fig. 7. 

Again, to observe the effect of mutual liabilities we perform an equivalent 
computations, but with zero mutual liabilities and the AP applied. The difference 
of two solutions is presented in Fig. 8. As one would expect, our results demonstrate 
a significant difference in the area close to X\ = fj,= , X 2 = ^ 2 , i.e. the effect of 
the mutual liabilities is pronounced in this area not only for the marginal survival 
probabilities, but for CDS prices as well. 


6 Pricing first-to-default (FTD) contracts 

The price of the FTD F\^ solves the following terminal boundary value problem 9 : 

Fi, t +CiF 1 =n, (34) 

Fi (t, 0, X 2 ) = 1 -Ri, Fi (t, Xi, 0) = 1 - R 2 . 

F\(t, Xi,X 2 t 00 ) = fi,oo(t,Xi), F 1 (t,X 1 too,X 2 ) = f 2 ,oo(t,X 2 ), 
F 1 (T,X 1 ,X 2 ) = 0ol (Xl ,x a )eD + ^l 1 (X 1 ,X 2 )£D 1 + P 2l(X 1 ,.Y 2 )GD 2 - 


For F 2 ,t this can be done by analogy. 
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Fig. 7: CDS prices Ci(X\, X 2 ) computed 
by using a Hundsdorfer-Verwer scheme. 



*i 


Fig. 8: The difference in CDS prices 
Ci(Xi, X 2 ) computed with and without 
mutual obligations using the FD method. 


Here /^ ?0 Q (t, JQ), i = 1,2 is the solution of the one-dimensional terminal bound¬ 
ary value problem 


dtfi,oo(t,Xi) + Cif itO0 (t,Xi) = (35) 

fi,oo(t,0) = 1 -Ri, fi,oo(t,oo) = -Si(T-t), 
fi,oo(T,Xi) = (1 - Ri)lxi<nf- 

As it could be seen, /i j00 (i, Xi) = ci t00 (t, Xi) given in Eq.(32). Also /3j = ft(Xi, X 2 ) 
in Eq.(34) is defined as 

Pi = l-m.m[R lT (l),R l ], (Xi,X 2 )eD u * = 1,2, 

Po = l-min[min[fl 2 ,T( 7 i).- R 2 ],min[^ 1 )T ( 72 ),iii]], (Xi,X 2 ) € D. 

Similar to the previous section, it can be shown that the solution of this problem 
reads 

nOO rOO 

Ftf^X'uXtimJ J ftl ( x lA)£[ D UD 2 jG(r , ,X 1 ,X 2 )dX 1 dX 2 (36) 

/ / 

PT rOO 1 pT /*oo 

-ft/ / gi(r -s,X 1 )dX 1 ds+-(l-R 2 ) ds G x , (t - s, Xi, 0)dXi 
Jo Jo ^ Jo Jo 

+ V l ~ Rl) j 0 ds J 0 G xAr -s,o,x 2 )dx 2 , 

where the Green’s function G (t',X i,X 2 ) again is given in Eq.(23). 


6.1 Numerical experiments 

For the same reason as before a direct implementation of Eq.(36) is impractical 
from the computational point of view. However, again a similar trick can be applied 
to significantly improve the accuracy in computation of the boundary integrals. 
We describe it in more detail in Appendix D. 
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Again we ran the same test as in the previous section with parameters of the 
model given in Table 1, and the same = 0.02. We used the same FD method 
to verify our solution, see the previous section for the description of the method. 
The FTD prices are shown in Fig. 9. 

In order to understand the effect of the mutual liabilities on FTD prices, we 
repeated this test, but with zero mutual liabilities and the AP applied. The dif¬ 
ference of two solutions is presented in Fig. 10. Same picture can be obtained by 
using our analytical approach. 

As in the previous cases mutual obligations significantly influence FTD prices 
especially in the area close to X\ = , A '2 = p-f. 




Fig. 9: FTD prices Fi(X\, X 2 ) computed 
by using a Hundsdorfer-Verwer scheme. 


Fig. 10: The difference in FTD prices 
F\ (X \, X 2 ) computed without and with 
mutual obligations using the FD method. 


We also present the difference in the CDS and FTD prices for the first bank 
computed with and without mutual obligations and maturity T = 5 years. These 
results are given in Fig. 11, 12. It is seen that with the increase of maturity the 
effects of the mutual obligations decreases and in the limit of very long maturities 
almost disappears. 

To illustrate how the terminal distribution of prices looks like in some partic¬ 
ular example (which was schematically given in Fig. 2) in Fig. 13 the difference 
Fi(T, Xi, X 2 ) — Ci(T, I 1 A 2 ) is presented as a function of (Xi,X 2 ). Obviously, 
-Fi(T, Xi, X 2 ) is positive in D 2 while C\(T, X±, X 2 ) vanishes there (the red box 
in the right bottom corner of the Figure). And they also differ in a part of the 
domain D (the left bottom corner). In other points of the computational domain 
the values of F[ (T, X \, X 2 ) and C\ (T, X-]. X 2 ) coincide with each other. 


7 Calibration 

The model described in section 3 has three unknown parameters: cri, 02 , p- We use 
CDS prices for two assets and the FTD price for both assets to calibrate these 
parameters. The calibration is done in Matlab using a simple non-linear least 
square approach where every given point (quote) is taken with the same weight. 
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Fig. 11: The difference in CDS prices 
C 1 (X 1 ,X 2 ) computed with and without 
mutual obligations, T=5 years. 


Fig. 12: The difference in FTD prices 
Ci(Xi,X 2 ) computed without and with 
mutual obligations, T=5 years. 



Fig. 13: The difference F\(T,Xi,X 2 ) — Ci(T, X±, X 2 ), T=1 year. 


In a test experiment we use all the parameters as in Table 1 and also use Ayo = 
300, A 2y o = 300, ci = ?2 = 0.05. Then we calibrate ui,u 2 ,P in the following way. 
First, we set ai = 0.3, a 2 = 0.4, p = 0.5 and compute the prices of CDS and FTD 
using our algorithm. This gives us the quotes Ci = 0.05, C 2 = 0.0583, F\ = 0.0583. 
Then we run the calibrator to make sure it converges to the same values oi a±, a 2 , p 
to validate self-consistence of our approach. 

Finally, since we investigate how strong the effect of mutual obligations on 
the parameters of the model is, in the second test we ignore mutual obligations 
and apply our AP as was discussed in the previous sections. The results of such a 
calibration are presented in Table 2. 

A typical time necessary to get the values of the parameters in Matlab is about 
10 secs with the objective function tolerance set to 10~ 4 . The corresponding time 
if the FD algorithm is used with the grid 70 x 70 points in space and time step 0.03 
is about 12 times slower. Certainly, for longer maturities this difference increases. 
The results for T = 5 years are also presented in Table 2. Here the computed quotes 
are C\ = 0.2579, C 2 = 0.3182, F± = 0.336. As it can be seen, accounting for the 
mutual obligations significantly affects the values of the calibrated parameters. 
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Table 2: Results of calibration 



T = 1 yr 

T = 5 yrs 1 

Test 

<71 

(72 

P 

(71 

(72 

P 

MO 

0.300 

0.400 

0.500 

0.300 

0.400 

0.500 

NMO 

0.2819 

0.4421 

0.4936 

0.3189 

0.4234 

0.2942 

Dif, % 

6.0372 

-10.5373 

1.2801 

-6.3108 

-5.8616 

41.1670 


8 Conclusion 

In this paper we consider interlinkage (mutual obligations) of banks and their 
influence on marginal survival probabilities as well as CDS and FTD prices of the 
corresponding names. We use a simple model where banks’ assets are driven by the 
correlated Brownian motions with drift. The choice of the model is dictated by the 
advantage to get all the results in a closed form, at least in the 2D case. A more 
sophisticated model withs assets driven by a general correlated Levy processes 
has been already considered in [15]. However, the present description is more 
transparent and allows one to better understand the nature of the effect, and also 
adds CDS and FTD prices to the picture. In the 2D case we also calibrated this 
model to some artificial market quotes and showed that the mutual obligations 
must be taken into account to get the correct values of the model parameters as 
they significantly influence the results of calibration. To the best of our knowledge 
these results are new. 

Another less important, but perhaps interesting result is a closed form solu¬ 
tion for the marginal survival probabilities for two assets driven by the correlated 
Brownian motions with drift. This solution yet was not given in the literature, so 
we present it in this paper. To understand a financial meaning of this solution, we 
underline that for big banks due to various regulators requirements their assets 
cannot drop down much below their liabilities, which means that their recovery 
rates R should be almost 1 (or, possibly, even exceed 1). In this case when com¬ 
puting, e.g., marginal survival probabilities, the domain D 12 in Fig. 2 becomes a 
positive octant. 
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A Solution of Eq.(7) 


We need to prove that Eq.(7) has a unique solution. The below discussion is an alternative to 
the solution of this problem given by [9], 

First, consider two extreme cases. If no banks default, then a;+7 jlji > 1, Vi S [1, N], 
Obviously, the solution of Eq.(7) is 7 j = 1. Vi S [1,1V]. 

If all banks default, then Eq.(7) transforms to the form 

F( 7 ) = 7 , 

where ^( 7 ) denotes the lhs of Eq.(7). This is a fixed point problem 10 that can be solved by 
the fixed point iterations method. A sufficient condition for local linear convergence of fixed 
point iterations is that the Jacobian J(F( 7 )) has to obey the condition 

|J(F( 7 ))I < I- (37) 

To prove this in our case, represent the Jacobian in the explicit form 



0 

hi 

hi 

■ ■ In 1 


0 

L21 

L31 . . 

■ l n1 

J(F( 7 )) = 

ll2 

0 

132 ■ ■ 

■ ■ In 2 

= (rW 

L 12 

0 

L32 • ■ 

■ L N 2 


llN 

h n 

h n ■ ■ 

.. 0 

\k= 1 / 

L\n 

L2N 

L3N ■ ■ 

. 0 


By definition for any matrix \M\ = \\rriij \|, i,j E [1,7V] 

N 

det(\M\) = E sgn(x) n m i>x ., 

xes N i =1 


where the sum is computed over all permutations \ °f the set Sn = 
all Lij > 0 we have 


N 


r n 


-1 


J(F( 7 )) < 


e n l '.*< 




XSSjv i=l 


Lfc=i 


[1, 2,..., N], see [4], Since 


(38) 


Now observe that the numerator in Eq.(38) is a sum of the products of the N elements, 
and each such a product i) is positive, and ii) has its vis-a-vis in the denominator. However, the 
denominator contains also some additional positive terms, for instance ^i-> anc ^ therefore, 

J(F( 7 )) < 1 - 

As an example, when N = 2 


\J(Fh))\ 


L21 L 12 

L\ + L 12 1/2 + 1/21 


< 1. 


Thus, we proved that the condition Eq.(37) is always satisfied. Therefore, by the Banach 
fixed-point theorem ([11] the map F( 7 ) —> 7 is a contraction mapping on 7 , and this implies 


10 This actually is a linear system of equations. However, we want to solve it using a fixed- 
point iterations method to later apply this technique to the general Eq.(7). 
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the existence and uniqueness of the fixed point since a unit cube where 7 is defined is a compact 
metric space. 

These two extreme cases naturally give rise to the idea of how to solve Eq.(7) in general 
by using a fixed-point iteration method. Given the vector 7 from the previous iteration, we 
check the condition a\ + yqgG?. < 1 for all i E [1,1V]. If for some i = k this condition 
is not satisfied, we put 7 *. = 1 and exclude the equation for 7 *. from Eq.(7). Otherwise, this 
equation remains in the system. After this step is completed, and for instance, M out of N 
variables 7 were set to 1, we solve Eq.(7) for the remaining N — M variables. The uniqueness 
of the solution and convergence of the fixed-point iterations follow from the above proof. 


B Closed form representation of the integral 

Due to various regulators requirements assets of large banks cannot drop below their liabilities, 
which means that their recovery rates R should be close to 1. In this case when computing, 
e.g., marginal survival probabilities, the domain D 12 in Fig. 2 becomes a positive quadrant. 
Finding an analytical solution for survival probability in a positive quadrant with non-zero 
drift is a long standing problem, that to the best of the authors’ knowledge was not solved 
yet. A relevant literature includes [7, 18, 22, 26] and references therein. 

From a technical prospective we want to compute the integral 

n 00 /*oo 

Q 1 (t',X' 1 ,X' 2 ) = / dX 1 / dX 2 G(r',X U X 2 ) (39) 

Jo Jo 

where the corresponding 2D Green’s function is given by Eq.(23). The closed form solution 
for this integral is known when the drift £ in Eq.(39) vanishes, see [22] and references therein. 
However, if £ ^ 0 the closed form solution is not known yet. Here we derive this representation 
in the form of series of generalized and confluent hypergeometric functions. 

First, using polar coordinates R,(f) we rewrite Eq.(39) in the form 

2 00 pvj r 00 2 

Qi(t', R' ,</>') = -e K sm(u n <p') / sin[y n (l>)d<f> / Re’ y ^ R e~ aR I v (/3R)dR (40) 

Jo Jo 


where 

« = - — ~ ~ ( x ',0), P = R' /$, liP) = (®2 + p0i) sin <j> + p 02 cos <fi. a=T. 

Next, we use the Gegenbauer expansion of the complex exponential of two variables in 
terms of the ultra-spherical (Gegenbauer) polynomials, [2] 

oo 

&ixs = r{v) (!)■" J2 ik (’' + k ) J -'+kO) c Ux), ( 41 ) 

fc =0 

where C^(x) are the Gegenbauer polynomials ([1]), and the parameter v can be arbitrary 
chosen. It can also be seen as a Neumann series ([24]) of the exponential e lxt . By changing 
variables s = iS in Eq.(41) the latter can be transformed to 

/ S'\ ~ v °° 

e-S* = r(v) J E(- 1 )' C (^ + k)I v+k (S)C" k (x)- (42) 

Substitution of this representation with S = (3R and x = — 7 into Eq.(40) yields 

Qi{t',R',4>') = ^e K r(z/)(^) ^2 (-l) M ( l/ + ^)sin {v n <j>') 

^ ' n =1 /x=0 

rvj p 00 2 

• / sm^ n 4>)C^(-y(4>)/P)d<t> / R}~ v e~ aR I„ n {pR)I v+tJl (pR)dR. 

Jo Jo 


( 43 ) 



24 


Andrey Itkin, Alexander Lipton 


For the sake of simplicity, it does make sense to choose v = 1, and then use the identity 

([ 10 ]) 


[°°e - aR2 I Vn {PR)I„+MR)dR = 2 -^n-M- 1 a -( 1 'n+M+ 1 )/ 2 / 3 l 'n+M r [(! + ^ + v n)/‘2] 

Jo r(/i + i)r(y n + 1 ) 


(44) 


3^3 


Mn + jU + 1 + p + 2 + p + 1 

2 2 2 . 

/i. + 1 + 1 p + + 1 


where 3 ^ 3 ( 01 , a 2 , 0 ( 3 ; 61 ,& 2 , 63 ; z) is a generalized hypergeometric function ([3]). 

Further, observe that at v = 1 the Gegenbauer polynomials become the Chebyshev poly¬ 
nomials of the second kind which admit the representation ([ 1 ]) 

[n/2] 

Unix) = ]T i-l) k C™- k (2x) n - 2k , 
k =0 


where [a;] is the floor function. Therefore, the first integral in Eq.(43) assumes the form 
[M/2] 

— 2fc ^ (~i^—k 

k Jo 


[M/2] 

h = V 2' / “ 2fc (-l ) k C£~ k / sin(z/ n (/>)(0i sin</< + 0 2 cos <f>) II ~ 2k d<j>, 

k=0 Jo 


(45) 


with 


01 = 


02 + p 0 \ p 0 2 


and be the binomial coefficient. 


P 9 P ' 

The integral in the rhs of Eq.(45) can be taken in closed form and reads 

y rzu 

sm(v n <f)){ 6 i sin <j> + O 2 cos 0 ) l! ~ 2 k d 0 = 

0 


(46) 


io 2 2 k-fi -1 
uj 2 {fi — 2 k ) 2 — n 2 n 2 


ai [bi.Fi in) + b 2 Fi(-n)] + a 2 [bi.F 2 (n) + b 2 F 2 (-n)] 


1 

/ 7r n \ 

,fc+ 

i 

/ 7rn \ 

- 1 

-m) 

- 1 

-Ml 

2 

Vo; / 


2 

V w / 

1 

/ 7rn \ 

,fc+ 

1 

/ 7m \ 

- 1 
2 

f- P) 

V a; / 

- 1 
2 

f- P) 

V id / 


a i 


i&2 


9i — i0 2 


2k — n 


, a 2 = e l7rn , bi = (—l ) 4 1 a;(^ — 2k) + 7 rn, z = 1,2, 


where 2 i*i(a, 6 , c, x) is confluent hypergeometric function ([ 1 ]). 

Although in Eq.(46) the integral is represented as a function of a complex argument, it 
could be shown that it is real. For example, 

f s\n(v n <f))( 6 i sin </> + 0 2 cos <j))d<f) = 2 ^^2 2 [(“l) n ( 0 i sin(cj) + 0 2 cos (a;)) - 0 2 ] 

Jo id — 7V n 


y rzu 

sin(z/ n </))( 0 i sin </> + 0 2 cos <ft) 2 d(f) = 

0 


— 4a; 3 (0 2 + 0 2 ) + 2tt 2 O^n 2 lu 


27r 3 n 3 — 87 r nw 2 

+ (—l) n cj [ 7 r 2 n 2 ((0 2 — 0|) cos(2a;) — 20 i0 2 sin(2cj)) — (0 2 + 0 2 ) (n 2 n 2 — 4cj 2 )] 


etc. 
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C Computationally efficient representation of Eq.(33) 

First, let us mention that the problem given in Eq.(27) is defined in the semi-infinite domain 
X\ G [0,oo), X 2 E [0, 00 ). However, for practical purposes this infinite domain is always 
truncated by some reasonably large value Mi, i = 1,2. Thus, we consider Eq.(27) with X 2 E 
[0, M 2 ]. Strictly speaking, this truncation will change the Green’s function representation ([23]), 
however the error should be small when X' 2 —> 00 , or in other words it is within the truncation 
error of changing the upper boundary from 00 to M 2 . 

To exactly match the boundary conditions in Eq.(27), we replace C\(t, X^ 2 ) with a new 
function 


Ci(t,xx 2 ) = Ci(t, XX 2 ) - - (i - 

Function C'i (/,. Y. X 2 ) solves the following problem: 

Ci, t (t,X 1 ,X 2 ) + CC 1 (t,X 1 ,X 2 ) = S(t,X!,X 2 ), 

Ci(t,X u 0) = Ci(t,0,X 2 ) = Ci(t,X 1 ,M 2 ) = Ci(t, 00 , X 2 ) = 0, 


V(t,x 1 ) 


(47) 


x 2 , 


M 


Ci(T, Xi,X 2 ) — ail(x 1 ,X2)e[D U D 2 ] M 

The solution of this problem is given by the formula ([23]) 
roo r 00 

C l (t',X[,X' 2 )= dX 1 / dX 2 C\(r',X\,X 2 )G{t', Xi, X 2 ) 

Jo Jo 


(48) 

J (! - • 


(49) 


pt roo roo 

+ / ds dX 1 / dX 2 S(r' — s,Xi,X 2 )G(t' — s,Xi,X 2 ). 

Jo Jo Jo 


At X 2 —> 0 and X 2 —> 00 due to the boundary conditions for C\ ( t , Xi, X 2 ) the new function 
Ci(t, X\, X 2 ) vanishes. Also, according to the boundary conditions ci j00 (t, X\) —> —$i(T — t) 
at X± —> 00 as well as ci,o(i, Xi), and Ci(t, Xi, X 2 ). Therefore, in this limit, Ci(t,Xi, X 2 ) 
vanishes as well. Finally, at X\ = 0 we have ci )O o(^0) = ci 5 o(f,0) = C\(t, X±, X 2 ) = 1 — 
R\. Therefore, in this limit, C\(t, X\, X 2 ) = 0. Thus, function C\(t, X\, X 2 ) satisfies the 
homogeneous boundary conditions. 

Now, let us give an exact representation of S(t, X\, X 2 ). We need to apply the operator 
dt + C to both parts of Eq.(47) and take into account Eq.(27) for Ci(t, Xi, X 2 ). Omitting a 
tedious algebra we obtain 


S(t, Xi ,X 2 ) = J2 a i (*> - Y 1. X 2) 

i =1 

a 1 (t,X 1 ,X 2 ) = <5(Xi -fi<)di(t,X 2 ), di(t,X 2 ) = 
= di(t)X 2 + d 2 (t) 


9ci,p(t, Xi) | 
dX~i 


lA'i=Af 


X 2 - M 2 

m 2 


a 2 (t, Xi,X 2 ) = <5'(Xi -/*<)a 2 (t,Xi,X 2 ), d 2 (t,X 1 ,X 2 ) = [ci,o(t, Xi) - ci,o(t, £<)] 


(50) 


x 2 - m 2 

2 M 2 ’ 


o 3 (t,Xi,X 2 ) = 1 Yi> -<a 3 (f,.Yi,.Y 2 ), a 3 (f,Xi,X 2 ) = [^ 2 (i) (ci,o(t, X 4 ) - )) 

+ (X 2 — M 2 ) (ft -9 tCl ,o(f,Af)) + o(t,Xi)] EE X 2 6i(t,X!) + 6 2 (t,Xi), 

o 4 (t,Xi,X 2 ) = a 4 (t, Xi) = -ft + J— [$ 2 (f) (ci j00 (f,Xi) - ci 

M 2 l 

+ p(d t c 1 , 0 (t,X 1 ) - 9iC lj00 (f,X 1 )) j. 

Further, denote 

y rT roo roo 

ds dX! dX 2 ai{r' — s,X 1 ,X 2 )G(t' — s,X\,X 2 ). 

0 Jo Jo 
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By using Eq.(49) and Eq.(50) we obtain 

/*t / r oo /* oo 

ds dX 1/ dX 2 E(r'- s,X 1 ,X 2 )G{t'-s,X 1 ,X 2 ) = JZ Ji 

Jo Jo Jo i=1 

PT POO 

Ji= ds 01 (r'- s,Xi,X 2 )G(t' - s, fi< ,X 2 )dX 2 

Jo Jo 

= f ds[d 2 (r'- s)Yi(r'- s,ii<) + di(r'- s)Y 2 (t'- s,fL<)\ , 

Jo 

PT POO 

J 2 = ds di(T’- s,^,X 1 ,X 2 )Gx 1 (t'- s,ix^,X 2 )dX 2 

Jo Jo 

= f ds \d 2 (r' — s)Z\(t' — s, jl^) + d\(r' — s)Z 2 {r' — s, &<)] , 

Jo 

PT POO POO 

J 3 = ds dX 1/ dX 2 a 3 {r' -s,X 1 ,X 2 )G{t’ -s,X u X 2 ) 

Jo j iif Jo 

PT ■- POO POO -■ 

= J ds[J dX\b 2 (r' — s,Xi)Yi(t' — s,Xi) + J dXi&i(r' - s,Xi)y 2 (r' - s,Xi)J, 

y PT POO 

ds / dXia 3 (T'-s,Xi)yi(T'-s,Xi), 

0 v/o 

where 

/* OO POO 

Y 1 (t,X 1 )= / G(t,Xi,X 2 )dX 2 , y 2 (t,Xi)= / X 2 G(t,Xi,X 2 )dX 2 , 

Jo v/o 

p 00 p 00 

Zi(t,Xi)= / G Xl (t,Xi,X 2 )dX 2 , Z 2 (f,Xi)= / X 2 G Xl (t,Xi,X 2 )dX 2 . 

v/o </o 

Also we emphasize that a pretty similar approach can be used for computing marginal 
probabilities. 


D Computationally efficient representation of Eq.(36) 


By using a similar idea as in the previous Appendix we first truncate the infinite domain 
(Xi, X 2 ) E [0, 00 ) x [0, 00 ) to a finite domain (Xi, X 2 ) E [0, M\\ x [0, M 2 \ and introduce a new 
function 


F±(t, X,X 2 ) 
h(t,Xi,X 2 ) 


Fi(t,XX 2 )-h(t,Xi,X 2 ) 



/l,oo(t,Xi) 



(1 - R 2 ) 


1 x l 


+ (1 — J?i) (l — lxi) 


~ /2,oo(*i V 2 ) [l — IM 1 -X 1 ] • 


(51) 


Function F\ (f . X } X 2 ) solves the following problem: 

F 1 ' t (t,X 1 ,X 2 ) + CF 1 (t,X 1 ,X 2 ) = T(t,X 1 ,X 2 ), 

Fi(t,Xi,0) = F 1 (t,0,X 2 ) = F 1 (t,X 1 ,M 2 ) = Fi(t, M\,X 2 ) = 0, 
F\(T, X\, X 2 ) = Fi(T, Xi,X 2 ) - h{T,Xi,X 2 ), 

and fi,oo(T,Xi) = (1 - Ri)l Xi <^=- 


(52) 
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The solution of this problem is given by the formula ([23]) 

POO POO 

F 1 (t',X' 1 ,X' 2 )= dX 1 / dX 2 F 1 (T',X 1 ,X2)G(T',X 1 ,X 2 ) (53) 

Jo Jo 

PT POO POO 

+ / ds dX 1 / ^Tfr' - s,Xi,X2)G(t' - s,Xi,X 2 ). 

Jo Jo Jo 

In order to compute T(t, Xi, X 2 ) apply the operator dt + £ to both parts of Eq.(51) and 
take into account Eq.(34) for Fi(t, Xi, X 2 ). Omitting a tedious algebra we obtain 

T(t,X i,X 2 ) = -2 ft + I [-1 + Ri + f2,oc(f ,X 2 )] 5’ Xl (M, - Xj) (54) 

5 ' ( 0 ) 

- \M2(Ri - R 2 ) - 2/(1 -R 2 + /i,oo(t,-Yi))] 

M 2 

+ 5(Xi)6i(t, Xi, X 2 ) + <5(Mi - Xi)6 2 (t,Xi,X 2 ), 

where 6i(£,Xi,X 2 ) are some functions. We omit the explicit form of these functions since the 
integrals 


I. 



5(AT)G(t' - 
<5(Mi - Xi)G( 


s,X 1 ,X 2 )b 1 (t,X 1 ,X 2 )dX 1 =0, 


r' - s, X 1 ,X 2 )b 2 (t, X x ,X 2 ) = 0, 


due to the boundary conditions for the Green’s function. Therefore, the final representation 
for the boundary integral in Eq. (53) reads 


PT POO POO 

ds dX i/ dX 2 T(r' — s, Xi, X 2 )G(t' — s, Xi, X 2 ) = K± + K 2 + K$, 

Jo Jo Jo 

pt' poo 

A'i = -2ft / ds dXili(r'- s,Xi) 

Jo Jo 

1 C T poo 

A'l = - / ds l-f 2 ,oo{t,X 2 ) + Ri - 1] G Xi (t’ - s, M\, X 2 )dX 2 , 

J Jo Jo 


K Z = 


r 

Jo 


(Ri - R 2 )Z\(t' - s,0) + Rl + - -Z 2 (r' - s, 0) 

M 2 


ds. 


At numerical (discrete) realization all Ki , i — 1 — 3 vanish at the boundaries as well as the 
first integral in Eq.(53), and so does Fi(t, Xi, X —2). Therefore, by definition of Fi(t, X±, X —2) 
this preserves the correct boundary conditions for Fi(t, X\, X — 2). 



